Apparent finite-size effects in the dynamics of supercooled liquids 
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Molecular dynamics simulations are performed for a supercooled simple liquid with changing 
the system size from N = 108 to 10 4 to examine possible finite-size effects. Although almost no 
systematic deviation is detected in the static pair correlation functions, it is demonstrated that the 
structural a relaxation in a small system becomes considerably slower than that in larger systems 
for temperatures below T c at which the size of the cooperative particle motions becomes comparable 
to the unit cell length of the small system. The discrepancy increases with decreasing temperature. 

PACS numbers: 64.70.Pf, 66.10.Cb, 61.43.Fs 
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As liquids are cooled toward the glass transition tem- 
perature T g , a drastic slowing-down occurs in dynamical 
properties, such as the structural relaxation time, the dif- 
fusion constant, and the viscosity @J§|, while only small 
changes are detected in static properties. The goal of 
theoretical investigations on the glass transition is to un- 
derstand the universal mechanism which gives rise to the 
drastic slowing-down. To this end, a great number of 
molecular dynamics (MD) simulations has been carried 
out for supercooled liquids j|. Several large scale sim- 
ulations have been performed very recently and reviled 
that the dynamics in supercooled liquids are spatially 
heterogeneous ^-^|; rearrangements of particle config- 
urations in glassy states occur cooperatively involving 
many molecules. We have examined bond breakage pro- 
cesses among adjacent particle pairs and found that the 
spatial distribution of broken bonds in an appropriate 
time interval (~ r Q ~ O.lrf, where r a is the structural 
a relaxation time and Tb is the average bond life time) 
is very analogous to the critical fluctuation in Ising spin 
systems. The structure factor is excellently fitted to the 
Ornstein-Zernike form j|, and the correlation length £ 
thus obtained grows rapidly with decreasing tempera- 
ture. Furthermore, we demonstrated that £ is related to 
t q trough the dynamical scaling law, t q ~ £ z with z ~ 4 
in 2D and z ~ 2 in 3D. The heterogeneity structure in 
our bond breakage is essentially the same as that in local 
diffusivity p], which leads to a systematic violation of 
the Stokes-Einstcin low in supercooled states. 

To investigate long-time behavior of glassy materials 
by MD simulation, rather small systems typically com- 
posed of N — 10 2 — 10 3 particles have been used with 
the periodic boundary condition (PBC). Such small sys- 
tems have generally been considered to be large enough 
to avoid finite-size effects in the case of amorphas mate- 
rials in which no long-range order exists. In fact, static 
properties such as the radial distribution function g(r) or 
the static structure factor S(q) of glassy materials are not 
seriously affected by the system size as long as reasonably 
large systems (N > 10 2 ) are used. However, this is not 
always the case for dynamical properties. For example, it 



is known that the use of a small system with PBC gives 
an manifest effect in relatively short-time behavior of the 
density-time correlation function. There appears an ar- 
tifact in time scale of order L/c, where L is the size of 
the simulation cell and c is the sound velocity Jl(]|ll| ■ As 
we already mentioned, the dynamical correlation length 
£ in supercooled liquids grows rapidly with lowering the 
temperature. It is thus possible that some kinds of finite- 
size effects may appear in the dynamics of supercooled 
liquids when £ becomes comparable to L even if no such 
effect is detected in the static correlation functions. The 
main purpose of this paper is to examine carefully this 
point for a simple soft sphere mixture. 

Our model mixture is composed of two soft sphere com- 
ponents 1 and 2 having the size ratio 01/02 = 1/1.2 
and the mass ratio 7711/777,2 = 1/2 while ei = £2 = e. 
The units of length, time, and temperature are o~i, tq = 
(mio-f/e) 1 / 2 , and e/ks in this paper. Details of simula- 
tions are given in our earlier paper ^| . We presently per- 
formed MD simulations only in three-dimensional space 
with the systems composed of N = Ni + N2 = 108, 
10 3 , and 10 4 particles while the density p = N/L 3 = 0.8 
and the composition Ni/N = 0.5 are fixed. The corre- 
sponding system linear dimensions are L N=108 = 5.13, 
l jv=io 3 = 1Q 8; and ^a^io 4 = 23.2. Simulations were 
carried out at T = 0.772, 0.473, 0.352, 0.306 and 0.267 
with the time step At — 0.005. The PBC was used in all 
cases. At each temperature, the systems were carefully 
equilibrated in the canonical condition so that no appre- 
ciable aging effect takes place. Data are then taken in 
the microcanonical condition. 

We first calculate the partial static structure factor, 
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to investigate whether finite-size effects are detectable in 
static particle configurations. Here r° and r\ are the 
positions of the j-th and A:-th particles in the a and b 
components (a, & 6 1, 2) and (■ • •) indicates the ensemble 
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average over different initial configurations. The dimen- 
sionless wave number q is in units of er^ 1 . In Fig. 1, we 
plotted S n (q) for N = 108, 10 3 , and 10 4 at T = 0.473 
(a) and 0.267 (b). One can find that Su(q) for all cases 
excellently agrees with each other both in (a) and (b); no 
systematic size dependence can be detected among them. 
We examined also S12 and S22 and confirmed the same 
tendency. Our results indicate that finite-size effects are 
very small or almost negligible in static pair correlations 
for N > 10 2 , as that is generally believed. 

Let us next consider finite-size effects in dynamical 
properties. The structure relaxation in glassy materi- 
als can be measured by calculating the coherent or the 
incoherent intermediate scattering functions, F(q, t) or 
F s (q,t). The decay profiles of those two functions tend 
to coincide at the first peak wave number q rn in F(q, 0). 
This has been confirmed for the present soft-sphere mix- 
ture H and also for a Lennard- Jones binary mixture p2| . 
Since F s (q, t) can be more accurately determined via MD 
simulation, we here calculate the incoherent scattering 
function for the component 1, 

Fs (q, t) = ±- expfr • Ar) (t)] J , (2) 

where Arj(t) = rj(i) — rj(0) is the displacement vector. 
Although F s (q,t) decays monotonically in normal liquid 
states, it exhibits multi-step relaxations in highly super- 
cooled states. This is due to the fact that at lower tem- 
peratures the particle motions are highly jammed and 
thus trapped considerably in effective cages formed by 
their neighbors. We then defined the a relaxation time 
T a , which corresponds to a characteristic life time of the 
effective cage, by F s (q,T a ) — e _1 at q = 2n for sev- 
eral temperatures. Fig. 2 shows the decay profiles of 
F s (q = 2w,t) obtained for N = 108, 10 3 , and 10 4 at 
T = 0.473 and 0.267. At T = 0.473, we see that the two 
curves from N = 10 3 and 10 4 entirely coincide, and one 
from N — 108 is also close to them. The relaxation times 
thus obtained are t^= 108 ~ 3.5, r^ =1 ° 3 = r^ =1 ° 4 ~ 2.3. 
However, the situation is different at T = 0.267 where 
F s (q,t) exhibits two-step relaxations. The faster and 
the slower parts of the decay are called the fast-/? (ther- 
mal) and the a relaxations, respectively. We note that 
the decay profiles for the three systems differ signifi- 
cantly in the a regime (t > 10 2 ), whereas they agree 
well in the fast-/3 regime (t < 1). Here we determined 

T AT=108 „ n000j T N=ltf „ 6500; and T N=1^ „ 2000 

at T = 0.267. Fig. 3 shows the temperature dependence 
of t^= 108 , r^ =1 ° 3 and r^ =1 ° 4 . At the highest tempera- 
ture T = 0.772, = 108 , r^ =1 ° 3 and t£ =1 ° 4 are exactly 
equal. However, t^ =108 begins to deviate from the oth- 
ers around T = 0.473 at which f ~ 5 Q is comparable 
to 7 J Ar = 108 — 5.13. The deviation increases with further 
decreasing temperature, and t^ =10S becomes almost one 
order larger than t^ =w for T < 0.352. Furthermore, 
t at=io 3 ]-, e g ms i deviate from r^ =10 around T — 0.306, 



at which L^^ 10 ' < £ < L Ar=1 ° 4 . We thus suppose that 
the present finite-size effects are attribute to suppressions 
of cooperative particle motions due to insufficient system 
size. The structural relaxation time of smaller systems 
thus tend to show a stronger (super- Arrhenius) tempera- 
ture dependence as a result of the finite-size effects. Re- 
membering the fact that the static structure factors are 
almost identical among those three systems at all temper- 
atures, the origin of this effect may be purely kinetic, or 
higher order correlations in particle configurations may 
relevant to this. In their recent paper |l3|, Horbach et 
al. found similar finite-size effects in a model silica glass 
which is known as a typical strong glass former, in addi- 
tion to the present soft sphere mixture which is usually 
classified in fragile glass former. 

To understand what happens in microscopic scale, we 
next visualize individual particle motions in N = 10 4 sys- 
tem at T = 0.267. First, we pick up mobile particles by 
the condition |Ar"(i)| > I" in a time interval [to, to + 1], 
where t = 0.125r Q = 250, and is defined separately 
for the component a e 1, 2 such that the sum of Ar°(t) 2 
over the mobile particles covers 66% of the total sum 
J2j Ar^(i) 2 . Then we define clusters of the mobile par- 
ticles by connecting i and j if |rj(t) — tv,-(0)| < 0.3(crj+(jj) 
or |j*{(0) — Tj(t)\ < 0.3(<Tj + (Tj) similar to Donati et al. 
||. In Fig. 4, we show spatial distribution of the clusters 
having the size n > 5; those are all chain- like |9| and have 
large scale correlations. Although only 5% of the total 
particles are shown in Fig. 4, the sum of Atv,-^) 2 cov- 
ers approximately 40% of the total Yl%=i ^ r i(*) 2 ' This 
clearly indicates that the cooperative motions become 
dominant in glassy states. To investigate finite-size ef- 
fects in cooperative motions quantitatively, we here in- 
troduce the distribution function, 

N IN 

P(n) = ' ^(tfS(n - m) Y, ' Ar ^ ' W 

i=l I i=l 

where the sum runs over mobile particles only, is the 
size of the cluster in which the mobile particle i belongs, 
and thus 6(n — n,) is 1 if i is a member of the cluster hav- 
ing the size n and if not. The physical meaning of P(n) 
is as follows; clusters having the size n contribute P(n) to 
the total squared displacements of the mobile particles. 
In Fig. 5, we show P(n) for N = 108, 10 3 , and 10 4 at 
T = 0.267 at which £ ~ 40 obtained for N = 10 4 is even 
larger than L N=1 ° = 23.2. We found that the coopera- 
tive motions in N = 108 system are strongly suppressed. 
By comparing N = 10 3 and 10 4 , it is found that larger 
scale cooperative motions (n > 10) are considerably sup- 
pressed also in N = 10 3 system. Infinitely large clusters 
which percolate the system though the PBC have never 
been found in all cases. The characteristic cluster size 
n = J2iLi nP(n) thus obtained is 3.38, 6.11, and 7.73 for 
N = 108, 10 3 , and 10 4 , respectively. 

In the framework of conventional liquid theories [ pd| , 
changes in static particle configurations upon cooling 
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lead to a drastic slowing-down toward the grass tran- 
sition. The mode coupling theory (MCT) is 
the most successful self-consistent approach within this 
framework. MCT describes onset of glassy slowing down 
(or slow structural relaxations) in the density-time corre- 
lation functions. In the original MCT, however, a sharp 
Ergodic/Non-ergodic transition is predicted at a temper- 
ature To which is considerably above T g . Although such 
a tendency has been found in colloidal systems in which 
thermal activation processes are negligible, it has never 
been observed in real glassy materials. It is thus be- 
lieved that the MCT has some difficulties for describing 
the true dynamics of supercooled liquids apparently be- 
low To. The main problem is the fact that the origi- 
nal MCT do not take into account the hopping motions 
of particles, which must be cooperative and thus long 
ranged as is seen in recent MD simulations . Unfor- 
tunately the problem has not yet been overcome in fully 
self-consistent way because efforts for including thermal 
activations make the theory more or less ad hoc. It is a in- 
teresting fact that the behavior of structural relaxations 
in our smallest system, in which cooperative hopping mo- 
tions are highly suppressed, becomes somehow closer to 
the original MCT prediction. 

It is worth mentioning several experimental attempts 
to find the evidence of the dynamical heterogeneity in 
glassy materials. One of the most interesting and useful 
approaches are the recent experiments on glass-forming 
thin films. The thickness d dependence of film properties 
is the main interest in these studies [^6|-^9|. The moti- 
vations of those studies are quite similar to the present 
study; they aimed to control the size £ of the cooperative 
motions by changing d and found that the relaxation time 
and T g considerably depend on d. They considered that 
the cooperative particle motions in the direction normal 
to the film may be truncated near the interface, and this 
effect become dominant when £ > d. Thus the system 
size restriction can enhance the particle motions. Note 
that this seems to contradict to the finite-size effect in the 
present MD simulations in which the size restriction sup- 
presses the cooperative particle motions |20| . The mech- 
anism of this discrepancy is still an open question; wc 
naively speculate that the situation is much more com- 
plicated in polymer films than in MD simulations. The 
system size restriction occurs only in one direction nor- 
mal to the film and other two in-film directions are free 
in thin films, whereas all directions are equally restricted 
in the present simulations. We believe that investigating 
the microscopic relaxation mechanisms in glass-forming 
(both simple and polymeric) thin films are definitely im- 
portant. 

In summary, we have examined system-size effects in 
the dynamics of supercooled liquids by MD simulation. 
We found significant finite-size effects in the structural re- 
laxation at lower temperatures, whereas no such effect is 
detectable in static pair correlation functions. The coop- 
erative particle motions, which leads to the a relaxation 
in glassy states, are strongly suppressed in smaller sys- 



tems for temperatures lower than T c at which £ becomes 
comparable to the system size. The present finite-size 
effects are regarded as a natural consequence of the dy- 
namical heterogeneity appearing in supercooled liquids. 
We point out that finite-size effects are significant in the 
dynamics of highly supercooled liquids, and thus special 
attention must be payed to investigate true relaxation 
dynamics in computer simulations particularly in the a 
regime. 
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FIG. 3. Temperature dependence of r Q for N = 108 (open 
diamonds), 10 3 (filled diamonds), and 10 4 (open squares). 
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FIG. 1. The partial static structure factor Su(q) obtained 
at T = 0.473 (a) and 0.267 (b) for N = 108, 10 3 , and 10 4 
systems. 




FIG. 4. Spatial distribution of particle displacements hav- 
ing the cluster size n > 5 at T = 0.267 for N = 10 4 . The 
arrows indicates individual particle displacements. 
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FIG. 2. The incoherent intermediate scattering function 
F s (q,t) of the component 1 with q — 2-n at T — 0.473 and 
0.267. 
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FIG. 5. P(n) versus n at T = 0.267 for N = 108 (open 
diamonds), 10 3 (filled diamonds), and 10 4 (open squares). 
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